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ESSA  RESEARCH  LABORATORIES 


The  mission  of  the  Research  Laboratories  is  to  study  the  oceans,  inland  waters,  the 
■WHflfc/ITttjfcu lower  and  upper  atmosphere,  the  space  environment,  and  the  earth,  in  search  of  the  under¬ 
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standing  needed  to  provide  more  useful  services  in  improving  man’s  prospects  for  survival 
as  influenced  by  the  physical  environment.  Laboratories  contributing  to  these  studies  are: 

Karth  Sciences  Laboratories:  Geomagnetism,  seismology,  geodesy,  and  related  earth 
sciences;  earthquake  processes,  interna!  structure  and  accurate  figure  of  the  Earth,  and 
distribution  of  the  Earth’s  mass. 


Atlantic  Oceanographic  Laboratories  and  Pacific  Oceanographic  Laboratories:  Ocean¬ 
ography,  with  emphasis  on  ocean  basins  and  borders,  and  oceanic  processes;  sea-air  inter¬ 
actions;  and  land-sea  interactions.  (Miami,  Florida) 


Atmospheric  Physics  and  Chemistry  Laboratory:  Cloud  physics  and  precipitation; 
chemical  composition  and  nucleating  substances  in  the  lower  atmosphere;  and  laboratory 
and  field  experiments  toward  developing  feasible  methods  of  weather  modification. 

Air  Resources  Laboratories;  Diffusion,  transport,  and  dissipation  of  atmospheric  con¬ 
taminants;  development  of  methods  for  prediction  and  control  of  atmospheric  pollu¬ 
tion.  (Silver  Spring,  Maryland) 


Geophysical  Fluid  Dynamics  Laboratory:  Dynamics  and  physics  of  geophysical  fluid 
systems;  development  of  a  theoretical  basis,  through  mathematical  modeling  and  computer 
simulation,  for  (he  behavior  and  properties  of  the  atmosphere  and  the  oceans.  (Princeton, 
New  Jersey) 

National  Hurricane  Research  Laboratory:  Hurricanes  and  other  tropical  weather  phe¬ 
nomena  by  observational,  analytical,  and  theoretical  means;  hurricane  modification  experi¬ 
ments  to  improve  understanding  of  tropical  storms  and  prediction  of  their  movement  and 
severity.  (Miami.  Florida) 

National  Severe  Storms  Laboratory:  Tornadoes,  squall  lines,  thunderstorms,  and  other 
severe  local  convective  phenomena  toward  achieving  improved  methods  of  forecasting,  de¬ 
tecting,  and  providing  advance  warnings.  (Norman,  Oklahoma) 

Space  Disturbances  Laboratory:  Nature,  behavior,  and  mechanisms  or  space  disturb¬ 
ances;  development  and  use  of  techniques  for  continuous  monitoring  and  early  detection 
and  reporting  of  important  disturbances. 

Aeronomy  Laboratory:  Theoretical,  laboratory,  rocket,  and  satellite  studies  of  the  phys¬ 
ical  and  chemical  processes  controlling  the  ionosphere  and  exosphere  of  the  earth 
and  other  planets. 

Wave  Propagation  Laboratory:  Development  of  new  methods  for  remote  sensing  of  the 
geophysical  environment:  special  emphasis  on  propagation  of  sound  waves,  and  electromag¬ 
netic  waves  at  millimeter,  infrared,  and  optica!  frequencies. 

Institute  for  Telecommunication  Sciences:  Central  federal  agency  for  research  and  serv¬ 
ices  in  propagation  of  radio  waves,  radio  properties  of  the  earth  and  its  atmosphere,  nature 
of  radio  no>so  and  interference,  information  transmission  and  antennas,  and  methods  for 
the  more  effective  use  of  the  radio  spectrum  for  telecommunications. 


Research  Flight  Facility:  Outfits  and  operates  aircraft  specially  Instrumented  for  re¬ 
search;  and  meets  needs  of  ESSA  and  other  groups  for  environmental  measurements  for 
aircraft.  (Miami,  Florida) 
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A  METHOD  OF  SMOOTH  CURVE  FITTING 
Hiroshi  Akima 


I 

A  new  mathematical  method  of  fitting  a  smooth 
curve  to  a  set  of  given  points  in  a  plane  is  developed,  and 
a  computer  subroutine  is  programmed  to  implement  the 
method.  This  method  is  devised  in  such  a  way  that  the 
resultant  curve  will  pass  through  all  the  given  points  and 
will  look  smooth  and  natural.  The  interpolation  between 
the  given  points  is  performed  locally,  and  no  assumption 
of  the  functional  form  is  made  for  the  whole  curve.  (  ) 

KT 

Key  words:  Direction  of  tangent,  interpolation,  poly¬ 
nomial,  smooth  curve  fitting.  4 


1.  INTRODUCTION 

When  we  try  to  determine  a  relation  between  two  variables,  we 
either  perform  computations  or  make  measurements.  The  result  is 
given  as  a  set  of  discrete  points  in  a  plane.  Knowing  that  the  relation 
can  be  represented  by  a  smooth  curve,  we  next  try  to  fit  a  smooth  curve 
to  the  set  of  points  so  that  the  resultant  curve  will  pass  through  all  the 
given  points.  Manual  drawing  is  the  most  primitive  method  for  this 
purpose  and  results  in  a  reasonable  curve  if  it  is  done  by  a  well-trained 
scientist  or  engineer.  But,  since  this  method  is  very  tedious  and  time 
consuming,  we  wish  to  let  a  computer  draw  a  curve.  To  do  so,  we  mu  it 
provide  the  computer  with  necessary  instructions  for  mathematically 
interpolating  points  between  the  given  points. 

There  are  several  mathematical  methods  of  interpolating  the 
value  of  a  function  from  a  given  set  of  values  (Milne,  1949;  Hildebrand, 
1956;  Ralston  and  Wilf,  1967),  but  the  application  of  one  of  these  methods 
to  curve  fitting  sometimes  results  in  a  curve  that  is  very  different  from 


one  drawn  manually.  In  other  words,  the  resultant  curve  sometimes 
appears  strange  and  unnatural.  In  this  report,  we  present  a  new  method 
of  interpolation  and  smooth  curve  fitting  that  is  devised  so  that  the  re¬ 
sultant  curve  will  look  smooth  and  natural. 

2.  DISCUSSION  OF  SOME  EXISTING  METHODS 

A  simple  example,  taken  from  a  study  of  FM  distortion  being 
conducted  by  the  author,  will  serve  to  illustrate  difficulties  encountered 
by  existing  mathematical  methods  of  interpolation.  Assume  that  the 
values  of  x  and  y  at  11  points  are  as  follows: 
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Knowing  from  the  physical  nature  of  the  ph<  omena  that  y(x)  is  a  single - 
valued  smooth  function  of  x,  we  try  to  interpolate  the  value  of  y(x)  and 
to  fit  a  smooth  curve  to  the  given  set  of  points. 

First  we  apply  the  method  of  interpolation  based  on  polynomials 
(Milne,  1949;  Hildebrand,  1956,  ch.  2,  3,  4).  This  method  is  perhaps 
the  one  most  often  used  because,  as  stated  by  Milne  (1949),  "polynomi¬ 
als  are  simple  in  form,  can  be  calculated  by  elementary  operations,  are 
free  from  singular  points,  are  unrestricted  as  to  range  of  values,  may 
be  differentiated  or  integrated  without  difficulty,  and  the  coefficients  to 
be  determined  enter  linearly."  There  are  several  variants  of  this 
method,  known  by  the  names  of  Newton-Cotes,  Lagrange,  Aitk/en,  and 
Neville.  Each  has  its  own  advantages  and  disadvantages, /but  they  are 
all  based  on  the  common  assumption  that  y(x)  can  be  closely  approxi¬ 
mated  by  a  polynomial  of  x  of  order  n  -  1,  where  n  is  the  number  of  given 
points.  They  should  give  the  same  result,  because  the  uniqueness  of  the 
polynomial  of  order  n  -  1  that  agrees  with  given  values  of  y(x)  at  the 
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given  n  points  has  been  proved  (Hildebrand,  1956,  p.  44).  The  result 
obtained  by  applying  the  10th  order  (Lagrangian)  polynomial  is  shown  in 
figure  1A  where  the  given  data  points  are  encircled. 

Next  we  try  to  apply  the  method  based  on  a  ratio  of  two  polyno¬ 
mials  or  a  rational  function  (Hildebrand,  1956,  sec.  9.9-9.12).  This 
method  is  not  so  commonly  used  as  the  polynomial  method.  Perhaps 
its  most  serious  disadvantage  is  that  the  desired  function  does  not  al¬ 
ways  exist;  in  our  example,  it  does  not.  Another  disadvantage  is  that 
the  non  singularity  of  the  function  is  not  guaranteed.  If  we  omit  the 
first  point  (0,  10),  however,  the  function  exists.  The  result  thus  ob¬ 
tained  is  shown  in  figure  IB. 

The  third  one  is  a  well-known  method  based  on  the  Fourier 
series  (Hildebrand,  1956,  sec.  9.  3),  which  exemplifies  the  so-called 
orthogonal  functions.  In  applying  this  method  to  our  data,  we  assume 
that  the  whole  range  of  x  from  0  to  10  corresponds  to  one-half  of  the 
fundamental  period  from  0  to  TT  and  apply  a  series  of  cosine  functions 
up  to  the  10th  order  harmonic  terms.  The  result  is  shown  in  figure 
1C. 

Finally,  we  apply  the  method  based  on  a  spline  function  (Ralston 
and  Wilf,  1967).  The  spline  function  of  degree  m  is  a  piecewise  func¬ 
tion  composed  of  a  set  of  polynomials,  each  of  order  at  most  m  and 
applicable  to  successive  intervals  of  the  given  data  points.  It  has  the 
characteristic  that  the  function  and  its  derivatives  of  order  1,  2,  . . . , 
m  -  1  are  continuous  in  the  whole  range  of  x.  This  function  includes 
the  (Lagrangian)  polynomial  as  a  special  case  when  m  =  n  -  1,  where 
n  is  the  number  of  given  points.  The  result  of  applying  the  spline  func¬ 
tion  of  degree  three  is  shown  in  figure  ID. 

In  addition  to  these  results  of  four  rather  representative  mathe¬ 
matical  methods,  the  curve  obtained  manually  is  shown  in  figure  IE. 
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It  is  the  average  of  curves  drawn  manually  by  several  scientists  and 
engineers.  A  comparison  of  the  curves  in  figures  1A  -  ID  with  IE  indi¬ 
cates  that  the  first  three  mathematical  methods  are  not  suitable  for  the 
example  given.  It  also  indicates  that  the  curve  obtained  by  the  spline 
function  and  shown  in  figure  ID  resembles  the  one  in  figure  IE,  but  we 
are  proposing  an  alternate  method,  discussed  in  the  sections  that  follow. 
As  we  see,  the  curve  obtained  by  the  new  method,  shown  in  figure  IF, 
is  closer  than  the  other  curves  to  the  manually  drawn  curve  in  figure 
IE. 

The  common  difficulty  encountered  by  the  existing  mathematical 
methods  is  that  the  resultant  curve  shows  unnatural  wiggles.  This  seems 
inevitable  if  we  make  any  assumption  concerning  the  functional  form  for 
the  whole  set  of  given  points  other  than  the  continuity  and  the  smoothness 
of  the  curve.  When  such  an  assumption  is  not  justified,  the  resultant 
curve  is  very  likely  to  behave  strangely,  as  in  figures  1A  -  ID.  The 
functional  form  is  what  we  try  to  derive  and  not  what  we  assume. 

When  we  try  to  fit  a  smooth  curve  manually,  we  do  not  assume 
any  functional  form  for  the  whole  curve.  We  draw  a  portion  of  the  curve 
based  on  a  relatively  small  number  of  points,  without  taking  into  account 
the  whole  set  of  points.  This  local  procedure  is  a  very  important  feature 
of  manual  curve  fitting  and  the  basis  for  our  new  method.  Note  that, 
although  the  spline  function  is  a  piecewise  function  composed  of  a  set  of 
polynomials,  all  polynomials  are  determined  simultaneously  on  the  basis 
of  the  assumption  of  continuities  of  the  function  and  its  derivatives  in  the 
whole  range,  and  no  individual  polynomial  can  be  determined  locally. 

It  is  not  easy  to  develop  a  mathematical  method  of  smooth  curve 
fitting  based  on  the  local  procedure.  If  one  of  the  existing  mathematical 
methods  is  applied  locally  or  piecewise  without  any  special  consideration, 
the  continuity  of  the  function  or  its  first  order  derivative  is  not  generally 
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guaranteed.  This  situation  is  illustrated  in  figure  2,  where  the  same 
data  as  in  figure  1  are  used  and,  for  each  N,  a  polynomial  of  order  N  -  1 
is  applied  to  each  set  of  N  successive  points  to  interpolate  the  part  of 
the  curve  in  a  unit  interval  around  the  center  of  the  N  points.  When  N 
is  an  odd  number,  the  resultant  curves  are  discontinuous;  when  it  is  an 
even  number,  the  slopes  of  the  resultant  curves  are  not  continuous,  and 
thus  the  curves  are  not  smooth. 

3.  NEW  METHOD 
3.  1.  Outline  of  the  Method 

Our  method  is  devised  so  as  to  work  in  two  different  ways,  i.  e. , 
one  for  a  single-valued  function  and  the  other  for  a  multiple -valued  func¬ 
tion,  to  correspond  to  the  two  ways  in  which  a  smooth  curve  is  fitted 
manually,  depending  on  whether  we  know  that  the  given  data  points  rep¬ 
resent  a  single-valued  function  or  not. 

Our  method  is  based  on  a  piecewise  function  given  by  a  polyno¬ 
mial  of  the  third  order  in  each  interval  for  the  case  of  a  single-valued 
function,  and  by  a  pair  of  polynomials  of  the  third  order  for  the  case  of 
a  multiple -valued  function.  For  a  single -valued  function,  our  method 
is  somewhat  similar  to  the  spline  function  of  degree  three.  As  in  the 
spline  function,  the  continuity  of  the  function  itself  and  of  its  first-order 
derivative  (the  direction  of  the  tangent  to  the  curve  or  the  slope  of  the 
curve)  are  assumed.  But,  instead  of  assuming  the  continuity  of  the 
second-order  derivative  as  in  the  third-degree  spline  function,  we  deter¬ 
mine  the  direction  of  the  tangent  locally  under  certain  assumptions.  By 
doing  so,  we  can  fit  a  curve  piecewise  to  the  given  set  of  data  points 
without  having  discontinuities  in  the  curve  and  its  slope. 

We  assume  that  the  direction  of  the  tangent  to  the  curve  (or  the 
slope  of  the  curve)  at  a  given  point  P»  is  determined  by  the  coordinates 
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of  five  points ;Pj_2 ,  P^,  Pt ,  P14.1,  and  P1+2.  In  other  words,  the  points 
more  than  two  intervals  away  are  assumed  not  to  affect  the  determination 
of  the  slope.  This  is  discussed  in  more  detail  in  the  next  section. 

The  portion  of  the  curve  between  a  pair  of  points  is  assumed  to  be 
determined  only  by  the  coordinates  of  and  the  slopes  at  the  two  points. 
This  interpolation  procedure  is  described  in  section  3.  3.  Since  the  slope 
of  the  curve  should  be  determined  also  at  the  end  points  of  the  curve, 
estimation  of  two  more  points  is  necessary  at  each  end  point.  This  ex¬ 
trapolation  procedure  is  described  in  section  3.4. 

3.  2.  Direction  of  the  Tangent 

Consider  five  pointSjl,  2,  3,  4,  and  5,  as  shown  in  figure  3A. 

Let  the  point  of  intersection  of  the  two  straight  lines  extended  from  line 
segments  12  and  34  be  denoted  by  A  and  the  same  point  corresponding 
to  line  segments  23  and  45  by  B.  We  seek  a  reasonable  condition  for 
determining  the  direction  of  the  tangent  CD  at  point  3. 

It  seems  appropriate  to  assume  that  the  direction  of  CD  should 
approach  that  of  23  when  the  direction  of  12  approaches  that  of  23,  and 
that  angle  /  23C  (the  angle  between  32  and  3C  )  should  be  equal  to  /  D34 
when  /  123  is  equal  to  /  345 .  With  these  rather  intuitive  reasonings  as 
a  guideline,  the  condition  of  determining  the  direction  of  CD  is  still  not 
unique.  For  simplicity  we  assume  that  the  tangent  CD  is  determined 
by  the  condition 

2~C  4D 

“  _  • 

CA  DB 

This  condition,  however,  does  not  exist  for  certain  configurations  of 
five  points,  such  as  the  one  shown  in  figure  3B.  In  this  case  the  alter¬ 
nate  condition, 
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2  C  _  4  D 

CA  DB 


does  exist,  as  shown  in  the  figure,  and  we  shall  use  this  condition. 

Summarizing  the  above  observations,  we  assume  that  the  direc¬ 
tion  of  the  tangent  CD  is  determined  by  the  condition 


CA  DB 


The  double  sign  depends  on  the  configuration  of  the  five  points,  and  the 
one  for  which  the  condition  exists  will  be  selected. 

Analytically,  as  shown  in  appendix  A,  condition  (1)  results  in  a 
quadratic  equation  with  respect  to  tan0,  where  9  is  the  angle  of  CD 
measured  from  the  x  axis.  For  a  given  configuration  of  the  points,  the 
discriminant  of  the  quadratic  equation  is  positive  or  zero  for  one  sign 
and  negative  for  the  other.  The  former  should  be  selected  as  a  matter 
of  course,  so  that  the  equation  has  real  roots.  Note  that  there  is  no 
ambiguity  in  the  selection  of  the  sign. 

When  one  line  segment  is  parallel  to  another  one,  the  direction 
of  the  tangent  CD  may  sometimes  be  indefinite  under  condition  (1). 
Analytically  this  corresponds  to  the  case  where  the  three  coefficients  of 
the  quadratic  equation  are  all  zero.  To  avoid  this  uncertainty,  we  as¬ 
sume  the  following: 

(a)  When  line  segment  23  is  parallel  to  34,  the  tangent  is 
parallel  to  23.  If  the  direction  of  23  is  opposite  that  of 
34,  the  tangent  changes  its  direction  by  180°  at  point  3. 

(b)  When  12  is  parallel  to  23  and  34  is  parallel  to  45,  the 
tangent  is  parallel  to  24. 

(c)  When  12  is  parallel  to  34  and  23  is  parallel  to  45,  the 
tangent  is  parallel  to  24. 
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Two  roots  or  solutions,  CD  and  C'D',  satisfy  condition  (1).  The 
one,  CD,  should  be  selected  so  that  the  two  points,  2  and  4,  lie  on  the 
same  side  of  the  straight  line  through  C  and  D.  This  can  be  done  by  add¬ 
ing  a  second  condition, 

sin  (/  C32)  •  sin  (/  C34)  >  0  .  (2) 

When  this  method  is  applied  to  a  multiple -valued  function,  the 
sense  of  the  direction  of  the  tangent  should  also  be  determined.  In  other 
words,  we  have  to  determine  both  cos  6  and  sin  9,  not  only  tan  0.  The 
sense  should  be  selected  so  that  both  the  positive  direction  of  the  tangent 
and  line  segment  34  lie  on  the  same  side  of  the  straight  line  through  2 
and  3.  This  can  be  done  by  adding  a  third  condition, 

sin  (0  -  G^j)  •  sin  (Ggy  -  G^y)  >0,  (3) 

where  G^j  and  G^t  are  the  angles  of  line  segments  23  and  34  measured 
from  the  x  axis,  respectively. 

Note  that  the  procedure  for  determining  the  slope  is  a  geometri¬ 
cal  one  and  is  independent  of  the  coordinate  system. 

3.  3.  Interpolation  Between  a  Pair  of  Points 
3.  3.  1.  Single -Valued  Function 

For  a  single -valued  function  y  =  y(x),  we  assume  that  the  curve* 
between  a  pair  of  points  can  be  expressed  by 

y  *  Po  +  Pi*  +  Pa*‘!  +  Pax3,  (4) 

where  the  p's  are  constants.  Since  the  coordinates  of  the  two  points, 
say  {xlf  yx )  and  (x2,  ys ),  as  well  as  the  directions  of  the  curve  tan  Gj 
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and  tan  6 s  at  the  points,  are  given,  we  further  assume  that  x  and  y  sat¬ 
isfy  the  conditions 


y  =  y-i  and  ^  =  tan  9j  at  x  =  xj 


y  =  y2  and  =  tan  Sa  at  x  =  xs , 


From  these  conditions  we  can  uniquely  determine  the  p  constants. 

Note  that  the  interpolated  curve  is  independent  of  the  scalings  of 
the  x  and  y  axes. 

3,  3.  2.  Multiple-Valued  Function 

For  a  multiple-valued  function,  we  assume  that  the  curve  be¬ 
tween  a  pair  of  points  (xj ,  yl )  and  (xa,  ys }  can  be  expressed  by 

x  =  Pc  +  Pi  z  +  Pbz"  +  poz3, 

(5> 

y  =  qo  +  qi  z  +  q2z2  t  q3 z3 . 

where  the  p's  and  q'»  are  constants  and  z  is  a  parameter  that  varies 
from  0  to  1  as  the  curve  is  traversed  from  (x; ,  yt )  to  (x8 ,  ya).  Since 
the  coordinates  of  the  two  points  (x; ,  y5 )  and  (x-  ,  ya  ),  as  well  as  the 
directions  of  the  curve  (cos  9  j. ,  sin  ©j, )  and  (cos  6a ,  sin  Sa )  at  the  points, 
are  given,  we  further  assume  that  x  and  y  satisfy  the  conditions 

dx 

x  =  Xi>  y =  yi-  di 

dx 

^  =  xa.  y  =  y=. 


=  r  r.os  9l(  and  =  r  sin©!  at  2  =  0, 

UZ 


=  r  cos  3a ,  and  =  r  sin  6a  at  z  =  1, 
dz 
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From  these  conditions  we  can  uniquely  determine  the  p  and  q  constants. 

Note  that  in  this  case  the  interpolated  curve  depends  on  the  seal' 
ings  of  the  x  and  y  axes, 

3.  4,  Extrapolation  of  the  Curve  at  an  End  Point 

Except  for  the  case  of  a  closed  curve,  estimation  of  two  more 
points  from  the  given  points  is  required  at  each  end  of  the  curve, 

3,4.  1,  Single -Valued  Function 

For  a  single -valued  function  y  =  y(x),  we  assume  that  the  curve 
near  the  end  can  be  expressed  by 

y  =  8o  +  gi*  +  gfcX3,  (6) 

where  the  g's  are  constants.  The  constants  can  be  determined  from 
the  coordinates  of  three  given  points  (xi ,  yx  J,  (xs,  ys),  and  (x^  ,  ya  ). 
Assuming  that 


X6  ■  X*  =  x*  -  Xj  =  Xa  -  xa  , 

we  can  determine  the  ordinates  y*  and  y6  corresponding  to  x*  and  xB , 
respectively,  from  (6), 

Note  that  the  extrapolated  points  are  independent  of  the  scalings 
of  the  x  and  y  axes. 

3.4.2.  Multiple -Valued  Function  (Noncloaed  Curve) 


For  a  multiple-valued  function,  we  assume  that  the  curve  near 


the  end  can  be  expressed  by 


x  =  Ho  +  Ki  z  +  . 

(7) 

y  =  ^  +  hji  t  hsz!  , 

when  the  g's  and  h'a  are  constants  and  z  is  a  parameter,  We  further 
assume  that 


x  =  Xj  and  y  =  y,  at  z  =  i  (i  =  1,  2,  3,  4,  5), 

From  the  coordinates  of  three  given  points  (x^ ,  yj ),  (xs ,  ys  ),  and 
(x0  ,  y3  ),  the  g  and  h  constants  can  be  determined,  and  consequently, 
also  the  coordinates  (x*  ,  y4  )  and  (xb  ,  yB). 

Note  that  the  extrapolated  points  are  independent  of  the  scalings 
of  the  x  and  y  axes . 


4.  EXAMPLES 

An  example  cx  the  application  of  this  method  is  shown  in  figure 
IF  where  the  curve  is  very  close  to  the  one  in  IE  determined  manually. 

Figure  4  gives  some  examples  of  the  application  of  the  method 
in  different  modes.  Two  curves  are  drawn  for  single-valued  functions 
y  c  y(x)  (MODE  =  1)  and  x  =  x(y)  (MODE  =  2)  in  A  and  B,  respectively. 
Two  examples  for  nonclosed-curve,  multiple -valued  function  (MODE  =  3) 
are  given  in  C  and  D,  A  circle  and  an  ellipse  are  drawn  in  E  and  F, 
respectively,  as  examples  for  the  case  of  closed  curve  (MODE  =  4). 

Figure  5  shows  artificial  examples  for  simple  configurations 
of  given  points  that  are  designed  to  supplement  the  description  of  the 
method,  especially  of  the  direction  of  the  tangent  (see  sec,  3.  2),  This 
figure  also  illustrates  how  strangely  curves  may  sometimes  behave  for 
adverse  configurations  of  given  data  points.  Application  of  the  third- 
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Figure  4.  Examples  of  application  of  the  new  method. 


Figure  5.  Further  examples  of  application  of  the  new  method. 


degree  spline  function  to  the  same  data  results  in  slightly  better,  but 
still  not  completely  satisfactory,  curves.  All  the  examples  in  this  fig¬ 
ure  are  periodic  functions,  but  neither  the  new  method  nor  the  spline 
function  method  is  devised  so  as  to  detect  or  take  advantage  of  regular¬ 
ity,  such  as  periodicity.  Because  periodicity  is  neglected  in  both  these 
methods,  the  resultant  curves  sometimes  appear  strange,  especially 
near  their  end  points. 

5.  CONCLUDING  REMARKS 

We  have  described  a  new  method  of  smooth  curve  fitting.  For 
proper  application  of  the  new  method,  thte  following  remarks  seem  per¬ 
tinent: 

(1)  The  curve  obtained  by  this  method  passes  through  all 
the  given  points.  Therefore,  the  method  is  applicable 
only  to  the  case  where  the  precise  values  of  the  coordi¬ 
nates  of  the  points  are  given.  It  should  be  recognized 
that  all  experimental  data  have  some  errors  in  them, 
and  unless  the  errors  are  negligible  it  is  more  appro¬ 
priate  to  smooth  the  data,  i.  e. ,  to  fit  a  curve  approxi¬ 
mating  the  data  appropriately,  than  to  fit  a  curve  pass¬ 
ing  through  all  the  points. 

(2)  Use  of  this  method  is  not  recommended  when  given 
data  points  manifest  apparent  regularity  or  when  we 
have  a  priori  knowledge  on  the  regularity  of  the  data. 

(3)  As  is  true  for  any  method  of  interpolation,  no  guaran¬ 
tee  can  be  given  of  the  accuracy  of  the  interpolation, 
unless  the  method  in  question  has  been  checked  in  ad¬ 
vance  against  precise  values  or  a  functional  form. 

(4)  The  method  yields  a  smooth  and  natural  curve  and  is 
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therefore  useful  in  cases  where  manual,  but  tedious, 
curve  fitting  will  do  in  principle. 

(5)  For  a  single-valued  function,  the  resultant  curve  is 
invariant  under  a  linear -scale  transformation  of  the 
coordinate  system.  In  other  words,  different  scal¬ 
ings  of  the  coordinates  result  in  a  same  curve. 

(6)  For  a  multiple -valued  function,  the  resultant  curve 
is  variant  under  a  linear-scale  transformation  of  the 
coordinate  system.  The  scalings  of  the  coordinates 
should  be  coincident  with  the  actual  size  of  the  graph. 

A  computer  subroutine,  named  CRVFIT,  has  been  programmed 
to  implement  the  method  reported  on  in  this  paper.  It  is  described  in 
detail  in  appendix  B. 
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APPENDIX  A 


Analytical  Expression  of  the  Condition  for  Determining 
the  Direction  of  the  Tangent 

Let  the  coordinates  of  points,  1,  2,  3,  4,  5,  A,  B,  C,  and  D  in 
figure  3  in  the  text  be  denoted  by  (xa ,  y!  ),  (x2 ,  y2),  (x3,  y3 ),  (x4,  y4 ), 
(xK,  ys),  (x»,  ya),  (xb,  yb ),  (xc ,  yc ),  and  (xd ,  yd ),  respectively,  and 
the  tangent  of  the  angle  of  CD  measured  from  the  x  axis  by  t.  Then, 
the  following  equations  should  hold: 


If  wc  introduce  new  constants  defined  by 

at  =  x1+1  -  Xj ,  (i  =  1,  2,  3,  4) ,  (14) 

bt  =  ym  -  yi,  (i  =  1*  2,  3,  4)  ,  (15) 

Sl3  ~  a,b4  -  ajbt,  (i  =  1,  2,  3;  j  =  2,  3,  4;  i<j),  (16) 
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and  eliminate  x's  and  y's  from  (8)  to  (13),  our  condition  reduces  to  a 
quadratic  equation  of  the  form 

At2  -  2Bt  +  C  =  0,  (17) 

where 

A  =  Sj2  Sg4  a3  T  S13  S34  ag  >  (18) 

B  =  Sjjg  Sg.-t  a3  b3  *F  Si3  S34  a2  b2  ,  (19) 

C  =  S12  S34  b3  T  S13  S34  b2  .  (20) 

The  discriminant  of  (17)  is  expressed  by 

D  =  B2  -  AC 

=  i  Ssl  S12  Sg4  Si3  Sg4  .  (21 ) 

We  will  take  the  upper  sign  in  (1)  and  (13)  when  the  product  Sls  S34  S13 
is  positive  or  zero,  and  the  lower  sign  when  it  is  negative.  By  doing  so, 
we  can  always  make  the  discriminant  nonnegative. 

If  the  coefficients  A,  B,  and  C  in  (17)  are  all  zero,  the  solution 
of  the  equation  is  indefinite.  This  occurs  when 


(a) 

S33  =  0, 

(b) 

s12  =  0 

and 

S34 

(c) 

S13  =  0 

and 

S34 

(d) 

s12  =  0 

and 

s13 

(e) 

S34  -  0 

and 

s34 

For  the  case  (a),  we  assume  the  provision  that  the  tangent  t  is  equal  to 
ba/a-2  =  b3/a3.  For  cases  (b)  and  (c),  we  assume  that  the  tangent  t  is 
equal  to  (b2  +  b3)/  (as  +  a3).  Cases  (d)  and  (e)  are  special  cases  of  (a), 
and  the  provision  for  (a)  applies. 
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I 


Equation  (17)  has  two  roots,  i.  e. , 

.  B  ±  JD 
1  "  A 


(22) 


Of  the  two,  the  desired  one  will  be  selected  so  that  it  satisfies 


Sao  Sqq  ^  0,  (23) 

where 

Sjjo  =  a2  bo  ~  ao  bg»  (24) 

Sco  —  ao  b3  ~  a3  b3,  (25) 

a0  =  1 N  1  +  tJ  ,  (26) 

bo  =  t//l  +  ty".  (27) 


Condition  (23)  can  equivalently  be  written  as  (2)  in  the  text. 
If  S20  satisfies 


S20  S33  >  0,  (28) 

cos  9  and  sin  9  are  equal  to  a0  and  b0,  respectively.  Otherwise,  cos  9 
and  sin  9  are  equal  to  -a0  and  -b0,  respectively.  Condition  (28)  is 
equivalent  to  (3)  in  the  text. 
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APPENDIX  B 


Computer  Subroutine  CRVFIT 

The  CRVFIT  subroutine  is  a  FORTRAN  subroutine  programmed 
for  the  CDC-3800  computer  to  implement  the  method  of  smooth  curve 
fitting  reported  on  in  the  text.  Necessary  information  for  the  user  of 
the  subroutine  is  given  on  pages  25  and  26  in  a  format  compatible  with 
the  library  function  manual.  A  FORTRAN  listing  of  the  subroutine  is 
given  on  pages  27  to  30.  A  binary  deck  of  the  subroutine  is  available 
from  the  library,  Computer  Division,  ESSA  Research  Laboratories, 
Boulder,  Colorado,  under  the  name  of  E2-IERB-CRVFIT. 

Although  the  CRVFIT  subroutine  has  been  programmed  for  the 
CDC-3800  computer,  it  can  be  modified  without  difficulty  for  other  digi¬ 
tal  computers  that  accept  a  FORTRA.N  language. 

The  CRVFIT  subroutine  is  devised  so  as  to  compute  coordinates 
of  closely  spaced  points  on  a  smooth  curve  determined  by  a  set  of  given 
points,  and  not  to  compute  a  value  of  ordinate  for  a  specific  value  of 
abscissa.  But  there  is  no  difficulty  in  modifying  the  subroutine  for  this 
purpose. 
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CRVFIT 


CRVFIT 


PURPOSE:  To  fit  a  smooth  curve  to  a  set  of  given  points  in  a  plane; 

i.  e. ,  to  compute  coordinates  of  a  new  set  of  points  (output 
points)  that  are  located  on  a  smooth  curve  determined  by 
a  set  of  given  points  (input  points)  in  a  plane  and  are  more 
closely  spaced  than  the  set  of  input  points  on  the  curve. 
(This  subroutine  interpolates  points  between  each  pair  of 
input  points,  and  the  output  points  consist  of  both  the  input 
points  and  the  interpolated  points.) 

FORTRAN  CALLING  SEQUENCE: 

CALL  CRVFIT(MD,L,X,Y,M,N,U,  V) 

where 

MD  =  mode  of  the  curve  (input  parameter), 

=  1  for  a  single-valued  function  Y  =  Y(X), 

=  2  for  a  single- valued  function  X  =  X(Y), 

=  3  for  a  multiple-valued  function,  nonclosed  curve, 
=  4  for  a  multiple-valued  function,  closed  curve, 

L  =  number  of  input  point  (input  parameter), 

X,Y  =  arrays  containing  the  abscissas  and  ordinates 
of  L  input  points  (input  parameters), 

M  =  number  of  divisions  between  each  pair  of  input 

points  (input  parameter), 

N  =  number  of  output  points  (output  parameter),  and 

U,  V  =  arrays  where  the  abscissas  and  ordinates  of 
N  output  points  are  to  be  displayed  (output 
parameters). 

ERROR  MESSAGE:  When  L  £  0  or  Ms  0,  the  error  message 

***  L  =  0/NEG  OR  M  =  0/NEG. 

A  =  (value  of  L)  Q  =  (value  of  M) 

ERROR  DETECTED  IN  ROUTINE  CRVFIT 

will  be  printed  on  the  standard  output  unit,  and  the  job  will  be 
aborted. 

STORAGE:  826  locations. 

TIMING:  (55M  +  710)  L  microseconds  for  MD  =  1, 

(68M  +  720)  L  microseconds  for  MD  =  2, 

(70M  +  1000)  L  microseconds  for  MD  =  3  or  4. 


(CRVFIT) 


METHOD: 

Step  1.  Except  for  MD  -  4,  two  more  points  are  estimated  at 
each  end  point  of  the  curve  by  assuming  that  the  curve 
near  the  end  point  can  be  expressed  by 

y  =  go  +  gi*  +  g2x2  for  MD  =  1 

and 

x  =  go  +  giz  +  g2Z2 
y  =  h0  +  ht  z  +  h2  z2  for  MD  =  3. 


Step  2.  The  direction  of  the  tangent 
to  the  curve  at  point  3  is 
determined  by  the  configu¬ 
ration  of  5  points,  1,  2,  3, 
4,  and  5,  shown  in  the 
figure,  by 


2  C 

4  D 

CA 

DB 

Step  3.  The  curve  between  a  pair  of  points  is  computed  by 
x  =  Po  +  Pi  z 

o  3 

y  =  qo  +  qi  Z  +  q2  z  +  q3  z  for  MD  =  1 

and 

x  =  po  +  pjz  +  p2z2  +  P3  z3 

2  3 

y  =  q0  +  qjZ  +  q2z  +  q3z  for  MD  =  3,  4. 

Computation  for  MD  =  2  is  made  by  interchanging  the  X  and  Y 
arrays,  applying  the  method  for  MD  =  1,  and  finally  interchanging 
the  U  and  V  arrays. 
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SUBROUT  INE  CKVF  1  f  t MOOE  *LO*X*Y.MO*NO*U*V» 
SMOOTH  CURVE  FITTING 

PROGRAMMED  oY  HIROSHI  AX  IMA*  ESSA-ITS 


MODE 
MODE 
MODE 
MODE 
LO  * 

wo  s 

MO  * 


*  l 

FOR 

Y  «  VIXI 

■  2 

FOR 

x  «  xtvr 

«  3 

FOR 

X  •  Xts>  AND 

v  » 

VIST  * 

NONCLOSED  CURVE 

■  4 

FOR 

X  *  A  IS)  AND 

Y  * 

VIST* 

CLOSED  CURVE 

NO. 

OF  INPUT  POINTS*  X 

AND 

Y  ARE 

INPOT  POINTS 

NO. 

OF  OUTPUT  POINTS*  U 

AND 

V  ARE 

OUTPUT  POINTS 

NO. 

OF  DIVISIONS 

C  UECLA3ATI0N  STATEMENTS 


DIMENSION 
DIMENSION 
tout VALENCE 

1 

EUUIVALENCE 

1 

2 

3 

4 

DIMENSION 
TYPE  DOUBLE 
t(JUl  VALENCE 

1 

DATA 


XI 101 *YI l0‘*UI loo f. VI  loo » 

A0»2»*80«2» 

( A  *AO 1 1 )  >  * IB  » AO« 2 ( *BO* 1  >  I  *  *  C *80 I  2  *  >  • 
IPO*X2>  *IQ0*Y2>  *«DX.A?> *IDV*B2  » 

I FLM*  TS*Z > • I JP* J5I • 

<DO.DA.D*Xl l • IDV*DB*R*ri t • 

IS2*S20*A1 »*IS3*S03*BI >* 

IP1.S12»*IP2*C12»*IP3*R12>* 

(Ql*S13>*(0?*C13!*:03*P13’ 

EWII2> *ER2»2) *MSGI3» 

DERI.DER2 

IL.E«1U)».IM.ERH2>  >.(DER1*ER1>« 
<LMI*ER2UH#«MM1*ER?I2»I*»DEP2*ER2> 
«MSG*2AHL  *  O/NEG  Oft  M  *  P/NEG.  I 


C  STATEMENT  FUNCTION 


SCRISIJ.CI J**ABSFtSiJ>-ABSFICIJS*1.0E-8 
C  PRELIMINARY  PROCESSING 


md=mode  * 

IFtL.LE.O.OR.M.LE.O) 

KP I =L  *M> 1  t 

DO  10  JP* 1 *L  * 

UOtPl )*XI IP>  ♦ 

10  CONTINUE  S 

DO  20  1*2 »L  % 

IFlUt<P2» .EO.UCKP3I.AND 
KP3  =  KP3-*-M  5 

20  CONTINUE  S 

IFIN.EO.l > 

IFtMD.NE.2) 

30  DO  40  KP4=1*N*M 

rs=u{KP4i  s 

40  CONTINUE 

SO  MM1*M-1  * 

IFIL.EO.2) 

LM1*L-1  S 


L*LPS*L0 
GO  TO  900 
IP*L-*1 
KP1=KPI-M 
VIKP1 >*YI IP> 
KP2*KP3*1 
KP2*XP2+M 
VIKP2I.EO.VIKP3M 
U«r.P3»*UIICP7) 
L*KP3/M*1 
GO  TO  890 
GO  TO  50 

U<KP4)*Vf  JCP4) 

FLM=M 
GO  TO  100 
GO  TO  200 


S 


$ 


* 

* 


% 

s 


M«M0 

IP* I P-1 


GO  TO  70 

VOCP3»*VOCP2) 

N*KP3 


V ( KP4 ) *TS 
DZ*1.0/FLM 
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C  SMOOTH  CURVE  FITTING  FOR  L  *  2 


IOO  DU*  tU(M) ~Ot 1 ) ) *DZ 

S 

DV*(VtN>-V(l> 1*02 

DO  Ho  K5*  X  *MMI 

UUS-*-l>*UUS)+DU 

ft 

VIKS+1  >*VIKS>*DV 

1X0  CONTINUE 

ft 

GO  TO  800 

SMOOTH  CURVE  FITT 

IMG  FOR 

L  GREATER  THAN  2 

200  X3*UU» 

ft 

V3*V*1 » 

X4-UIM+U 

ft 

74*VtM*l > 

K5*1*M+M 

ft 

X5-UU5* 

ft 

Y5*VCX5» 

A3-X4-X3 

ft 

B3*V4-V3 

A4*X5-X4 

ft 

B4*VS-Y4 

S34»A3*B4-A4*b3 

ft 

C34-A3*A4+B3*B4 

ft 

RS4*SCR(S34tC34) 

IFtMD#LE.3> 

GO  TO  400 

XI-N-M-M 

ft 

Xl-UtKl) 

ft 

Y1«VIX1) 

K.2-KI+M 

ft 

X2-UCK2) 

ft 

Y2*VU2> 

A1-X2-XI 

ft 

B1-V2-Y1 

A2-X3-X2 

ft 

B2-Y3-V2 

230  S12*AI*B2-A2»B1 

ft 

C12»A1*A2*B1*B2 

ft 

R12-SCRCS12.C12) 

S13*A1«B3-A3*B1 

ft 

Ci3»Al*A3*Bl*B3 

ft 

R13*SCR<Sl3tC13) 

S23"A2*B3— A3*B2 

ft 

C23*A2»A3+B2«B3 

ft 

R23=SCR(S23»C23) 

S24*A2*B4-A4»B2 

ft 

C24*A2*A4*B2»B4 

ft 

R24*SCRtS24,C24) 

ASSIGM  240  TO  LBL 

ft 

GO  TO  500 

240  DO  290  I-2.L 

ft 

X5*K5*M 

X2*X3 

ft 

Y2-V3 

X3*X4 

ft 

Y3*Y4 

X4-X5 

ft 

Y4*V5 

A2*A3 

ft 

B2*B3 

A3*A4 

ft 

B3*B4 

SI2*S23 

ft 

RI2-R23 

S13-S24 

ft 

RI3*R24 

S23*S34 

ft 

C23-C34 

ft 

R23*R34 

COS2-SGM»COS3 

ft 

SIN2*SGN*SIN3 

IF ( I.LT.LM1I 

GO  TO  270 

IF  I MD*EQ»4 1 

GO  TO  260 

IF  t I.EQ.LMI» 

GO  TO  450 

A4*A5 

ft 

B4*B5 

ft 

GO  TO  280 

260  IFIK5.GT.N) 

X5*l+M 

270  X5-UCK5! 

ft 

Y5*V(K5> 

A4-X5-X4 

ft 

94*Y5-Y4 

280  S24*A2»B4-A4*B2 

ft 

C24»A2*A4+B2*B4 

ft 

R24*SCR(S24  tC24» 

S34*A3*B4-A4»B3 

ft 

C34*A3*A4+B3*B4 

ft 

R34*SCR ( S34  »C34  J 

ASSIGM  600  TO  LBt- 

ft 

GO  TO  500 

290  CONTINUE 

ft 

GO  TO  800 
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EXTRAPOLATION  AT  THt  BEGINNING 


C 


400 

fFiR34.Lfc*0*G> 

GO  TO  430 

IF (MD.tQ.3) 

GO  TO  420 

4 1 Q 

Al*A2-A3 

S 

DBa2*A3*I»34/IA4»  ( A3.A4 )  ) 

b2=b3-D8 

* 

b 1 *B2-DB 

S 

GO 

TO 

230 

420 

UA= A4-A3 

s 

DB*B4-B3 

A2*A3-bA 

s 

B2-B3-DB 

A1*A2-L)A 

s 

8-1  *B2-DB 

% 

GO 

TO 

230 

43o 

A1*A2=A3 

s 

B 1 =B2*B3 

t 

GO 

TO 

230 

EXTRAPOLATION  AT 

THE  END 

450 

IF  t  R23-LE .0»0 ) 

GC  TO  480 

IF ( MD.t0.3 ) 

GO  TO  470 

460 

A5=A4»A3 

S 

D8  =  2«A3*S23/«A2*IA2+A3>  » 

b4*b3>Db 

* 

B5*B4.DB 

s 

GO 

TO 

280 

470 

UA«A3-A2 

% 

DB*B3-B2 

A4*A3>DA 

% 

Ba*B3*DB 

A5S  A4-.UA 

s 

B5=b4.DB 

s 

GO 

TO 

280 

480 

A5= A4=A3 

s 

B5»B4»B3 

% 

GO 

TO 

280 

OETERMlNAT IUN  OF 

THt  DIRECTION 

500 

5GN»1.0 

IFCR23.LE.O-0) 

GO 

TO 

550 

IF(R12.LE.0.0.AN0.R34. 

LE 

.0.01 

GO 

TO 

580 

IFCR13.LE-0.0-AND.R24. 

LE 

•  0.0) 

GO 

TO 

580 

IFCR12.LE.0.0.OR 

•  R24. 

LE 

•  0.0) 

GO 

TO 

560 

IFtR13.LE.0-0. OR 

.  R34. 

LE 

.0.0) 

GO 

TO 

570 

i>2  =  512*S24 

» 

63«S13*S34 

IF  (52#i>3.LT.0.0) 

S3*-S3 

A*S2»A3*A3-S3»A2 

*A2 

0=S2»A3*B3-S3*A2 

*B2 

C=S2*B3*B3-53»B2 

*B2 

D=S23*S0RTFIS2*S3 > 

IF(B*D.LT.0.0» 

D  =  -D 

» 

B*B+D 

S20  =  A2*B0«  H-A0« 

1  )*B2 

503*A0‘ 1)*B3-A3» 

80<  1  > 

IFCS20*S03.LE.0. 

0) 

GO 

TO 

510 

C0S3=A0I 1 > 

$ 

SIN3*B0( 1 » 

s 

GO 

TO 

520 

510 

S20=A2*B0«2>-A0( 

2)*B2 

COS3=AO<2) 

1 

SIN3*B0I2) 

520 

IF<S20*S23.GT.0. 

0> 

GO 

TO 

590 

COS3*-COS3 

* 

5 I N3=-S IN3 

s 

GO 

TO 

590 

550 

IF  tC23.LT .0*0  ) 

SGN=-1 .0 

560 

COb3=A 2 

$ 

SIN3-B2 

* 

GO 

JO 

590 

570 

COS3=A3 

s 

SIN3-B3 

* 

GO 

TO 

590 

580 

COS3=A2.A3 

s 

SIN3=B2.03 

590 

IFIMD.LE.2* 

GO  TO  LBL. 

(240.600) 

R  =  SUKTF(COS3#COS3-*-SIN3»SIN3) 
C053=COS3/R  *  5 1 N3=S IN3/R 

GO  TO  LBL  *  (240*600) 
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interpolation  in  a  section 

600  KSO=< I-2)*M+1  $  2=0.0 

I F I MD.GT . 2 )  GO  TO  660 

610  KS1 =KS0 
P1  =  0X 

Ql=PI*SIN2/COS2 

Q2= 3 *0Y-2*Ol-Pl*S I N3/COS3 

Q3=DY-01-02 

DO  620  JS=1*MM1 

KSl=KSl+l  $  Z=Z+DZ 

uiicsi  i=po+z*pi 

620  VOC.S1  )  =Qo+Z*  (Ol  +  Z*  <  Q2  +  Z*Q3  >  ) 

660  KS2=ICS0 

r=sortfiox*dx+dy*dy ) 

Pl=R*COS2 

P2= 3* DX-R* ( 2*COS2+COS3 ) 

P3=DX-P1-P2 
01 =R*S I N2 

Q2=3»DY-R*(2*SIN2+SIN3) 

03=DY-01-Q2 
DO  670  JS=1»MM1 

K.S2=KS2+1  $  Z  =  Z+DZ 

Ut<S2)=P0+Z*IPl+Z*(P2+Z*P3) ) 

670  VUS2)=Q0+Z*(Ul+Z*(Q2+Z*O3M 

NORMAL  return 

000  IFIMD.NE.2I 

DO  810  KR= 1 tN 
TS=U<KK> 

810  CONTINUE 
890  LO=LPS 

ERROR  EXIT 

900  DER2=DER1 
END 


GO  TO  890 

i  U<  KR ) =V( KR )  $ 

*  NO=N  $ 

S  CALL  08QERRORIO*MSG» 


GO  TO  290 


GO  TO  290 


V(K.R)=TS 

RETURN 
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